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We apply the Barker-Henderson (BH) perturbation theory to the study of a Lennard-Jones fluid confined in a 
random porous matrix formed by hard sphere particles. In order to describe the reference system needed in 
this perturbation scheme, the extension of the scaled particle theory (SPT) is used. The recent progress in the 
development of SPT approach for a hard sphere fluid in a hard sphere matrix allows us to obtain very accurate 
results for thermodynamic properties in such a system. Hence, we combine the BH perturbation theory with 
the SPT approach to derive expressions for the chemical potential and the pressure of a confined fluid. Using 
the obtained expressions, the liquid-vapour phase diagrams of a LJ fluid in HS matrix are built from the phase 
equilibrium conditions. Therefore, the effect of matrix porosity and a size of matrix particles is considered. 
It is shown that a decrease of matrix porosity lowers both the critical temperature and the critical density, 
while the phase diagram becomes narrower. An increase of a size of matrix particles leads to an increase of 
the critical temperature. From the comparison it is observed that the results obtained from the theory are in 
agreement with computer simulations. The approach proposed in the present study can be extended to the 
case of anisotropic fluid particles in HS matrices. 

Key words: fluids in random porous media, Barker-Henderson perturbation theory, liquid-vapour 
coexistence, scaled particle theory 
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1. Introduction 

The original understanding of the nature of the liquid state of matter is connected with Van der Waals 
equation of state formulated nearly 150 years ago (TJ. The Van der Waals picture focuses on different roles 
of the strong short-ranged repulsive and long-ranged attractive intermolecular interaction in forming the 
equilibrium properties of dense fluids. According to this picture, the harsh repulsive interactions fix the 
shape and size of molecules and essentially determine a high-density fluid structure. While the contribu¬ 
tion of repulsive interaction is entropic, the contribution of attractive interaction is mainly energetic and 
can be treated as the perturbation. The first theory of liquids based on the Van der Waals idea was pro¬ 
posed by Barker and Henderson (BH) nearly fifty years ago 12—1 j. Within this theory, the intermolecular 
potential is separated into the repulsive and attractive parts. The short-ranged repulsive contribution is 
described within the framework of the hard spheres model while the long-ranged attractive part is in¬ 
cluded within the Zwanzig high-temperature perturbation theory ( 3 ). A few years later, Andersen, Chan¬ 
dler and Weeks (ACW) developed a somewhat different theory of liquid based on the Van der Waals ap¬ 
proach (6}{9). In this theory, instead of intermolecular interactions, they separated intermolecular forces 
into repulsive and attractive parts. They also used the optimized cluster expansion (OCE) instead of the 
high-temperature perturbation theory for the treatment of attractive forces. However, the first term of 
perturbation related to the high temperature approximation (HTA) in both theories is identical except 
that the repulsive and attractive interactions are not exactly the same. 

‘Dedicated to Prof. Douglas Henderson on the occasion of his 80th birthday. 


©M.F. Holovko, T.M. Patsahan, V.I. Shmotolokha, 2015 


13607-1 






Holovko et al. 


In this paper we make use of the Van der Waals ideas and expand BH theory for the description of 
liquids adsorbed in random porous media. To this end, we use the Madden-Gland model (10). According 
to this model, a porous medium is presented as a quenched configuration of randomly distributed hard 
spheres forming a so-called matrix. A specific description of a fluid in such porous media is connected 
with double quenched-annealed averages: the annealed average is taken over all fluid configurations and 
the additional quenched average should be taken over all realizations of a matrix. One of the popular ap¬ 
proaches to the solution of this problem is based on the replica method, which allows for the extension 
of many theoretical methods of liquid state physics to the case of a fluid confined in a random porous 
medium. For instance, using the replica Ornstein-Zernike (ROZ) integral equation theory (TT), the statis¬ 
tical mechanics approach of liquid state was extended to the description of various fluids confined in 
random porous media fl2llT3l . including the chemically reacting fluids (Ml[IS). However, unlike bulk 
fluids, no analytical result has been obtained from the ROZ integral equation approach even for such a 
simple model as a hard sphere (HS) fluid in a HS matrix. At the same time, the model of HS fluid has 
a peculiar importance, since similar to the case of a bulk fluid COS 23. it can be used as the reference 
system in the development of different perturbation schemes. 

The first rather accurate analytical results for a HS fluid in a HS matrix were obtained quite recently 
in 1 11811231 by extending the scaled particle theory (SPT) 125|25) to the case of a HS fluid confined in 
random porous media. The SPT approach is based on the combination of the exact treatment of a point 
scaled particle in a HS fluid combined with the thermodynamical treatment of a macroscopic scaled par¬ 
ticle. The exact result for a point scaled particle in a HS fluid in a random porous medium was obtained 
in Cl. However, the approach proposed in Cl referred to as SPT1 contains a subtle inconsistency ap¬ 
pearing when the size of matrix particles is much larger than the size of fluid particles. Later on, this 
inconsistency was eliminated in a new approach referred to as SPT2 (20). Starting from this formalism, 
a series of new approximations were developed 12011231 . Among these approximations we select only 
SPT2b and SPT2bl, which will be used in this paper for the description of the reference system in the BH 
perturbation theory. 

The liquid-vapour phase diagrams play a central role in understanding the nature of liquids. In con¬ 
trast to the bulk case, the results of investigations of liquid-vapour phase equilibrium of a simple fluid in 
random porous media are rather controversial. Computer simulations of a simple fluid confined in a HS 
matrix I26il29l demonstrate a possibility of the existence of two phase transitions. One is analogous to the 
bulk liquid-vapor transition, but with a narrower coexistence curve and lower values of the critical den¬ 
sity and critical temperature. The second transition occurs at lower temperatures and at higher densities, 
and it is interpreted as a phenomenon related to the wetting effects in the fluid located in more confined 
regions of the matrix. On the other hand, in more thorough investigations (29), it was noticed that this 
second transition is extremely sensitive to a particular matrix configuration. In the case of HS matrix, 
it was shown that two phase transitions appear in some realizations of the matrix, while a single phase 
transition was observed in the others. In order to explain this observation, optimized cluster expansions 
were used in (30). It was found that different approximations can lead to qualitatively different results. 
For example, the mean spherical approximation (MSA) gives only a single liquid-vapour transition, but 
the inclusion of the second and the third cluster coefficients result in two phase transitions. At the same 
time, it was shown in (31) that the association theory leads to one phase transition. Our recent investiga¬ 
tion, in which we generalized the Van der Waals equation for simple fluids in random porous media, also 
demonstrates only one liquid-vapour transition (22). 

In the present study we combine the BH theory with the previously developed SPT approach to de¬ 
scribe a liquid-vapour phase behaviour of a Lennard-Jones fluid confined in a random matrix formed by 
hard spheres (HS matrix). For a comparison, the results obtained in different approximations are consid¬ 
ered as well. As it was mentioned above for the reference system, the SPT approach is applied using the 
SPT2b approximation and its improved version SPT2bl. It should be noted that the SPT2b approximation 
is our first really successful result for confined fluids (20), although it may have essential problems at 
high fluid densities and/or for low matrix porosities. At the same time, the improved SPT2bl approxima¬ 
tion is based on the original SPT2b, but it is free of this shortcoming and gives an accurate description 
of the thermodynamics for a hard sphere fluid in a hard sphere matrix up to the close packing condi¬ 
tions I2lti23l|32l . Apart from the BH theory, we present the results obtained in the HTA approximation. 
Using the developed approach, the phase diagrams for a confined Lennard-Jones are built. Different ma- 
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trix porosities and matrix particle sizes are considered. To check the accuracy of the approaches proposed 
in this paper, a comparison with the MSA results as well as with the results of Monte Carlo simulations is 
performed. 


2. Theory 


2.1. Reference system: HS fluid in HS matrix 


We start our theoretical consideration with the description of the reference system. For this purpose 
we briefly recapitulate the main ideas of the SPT theory and present here the expressions for the chemical 
potential and pressure of a HS fluid in a HS matrix, which were obtained in our previous papers I20II21I . 
and which are needed in the current study. The key point of the SPT theory consists in a derivation of the 
excess chemical potential of an additional scaled particle of a variable size inserted in a fluid. This excess 
chemical potential is equal to a work needed to create a cavity in a fluid which is free from any other 
particles. For a small scaled particle in a HS fluid in the presence of a porous medium, the expression for 
the excess chemical potential is equal to ( 20 ): 


/3ps X = lnp 0 (A s )-ln 


Cl + A s ) 3 
771 Pods) 


( 2 . 1 ) 


where p = l/(fc B T), k B is the Boltzmann constant, T is the temperature, 771 = ^npia^ is the fluid packing 
fraction, p\ is the fluid density, o\ is the diameter of HS fluid particles. The term po(A s ) = exp{-pp®) is 
defined by the excess chemical potential of the scaled particle confined in an empty matrix, p°. It has the 
meaning of probability to find a cavity created by the scaled particle in the matrix in the absence of fluid 
particles. We should note that here we use conventional notations l flTKl5]|T8l{22) . where the index “1” is 
used to denote a fluid component, the index “ 0 ” denotes matrix particles, while for the scaled particles 
the index “s” is used. 

For a large scaled particle, the excess chemical potential is presented by the thermodynamic expres¬ 
sion for the work needed to create a macroscopic cavity inside the fluid, which at the same time is con¬ 
fined in a porous medium, and the corresponding expression can be presented as follows: 

P 14 

Ppf=M A s ) + p—±~, (2.2) 

Po(A s ) 


where P is the pressure of fluid, V s is the volume of a scaled particle. The multiplier l/p 0 (A s ) appears 
due to an excluded volume occupied by matrix particles. In this context, it should be mentioned that the 
probability po(A) is directly related to two different types of porosity I20H22I I. The first one corresponds 
to the case of A s = 0 and provides the geometrical porosity 


00 = po(A s = 0), 


(2.3) 


which depends only on the structure of a matrix and it is equal to the volume fraction of a void between 
the matrix particles. For a HS fluid in a HS matrix, it is equal to 


0o = l-Po- (2.4) 

where rjo = ^na^po, p 0 = ^r, TVo is the number of matrix particles, ap is the diameter of the matrix 
particles, V is the volume of the system. 

The second type of porosity corresponds to the case A s = 1 and provides the probe particle poros¬ 
ity ( 2 QH 22 ) 

0 = p 0 (A s = 1) = e - ^°, (2.S) 

which is defined by the excess chemical potential of fluid particles in the limit of infinite dilution p®. 
Using the SPT theory |'24l[25) for the case of a HS fluid in a HS matrix, the following expression for 0 can 
be derived: 

0 =■ ( 1 - 770 ) exp 


3770T 3770(1 + ^770) -r pPpTjo 


1-770 


(1 -T7or 


P 0 


( 2 . 6 ) 
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where t - a 1 la 0 , and P 0 is the bulk pressure of the matrix particles. 

PPq = (1 + t?o + ??q) 

Po (l-??o) 3 

According to the ansatz of SPT I 18tf22lf2-h 251. w(A s ) can be presented in the form of an expansion: 

1 2 

w(A s ) = u>o+ u>iA s + - W 2 A s . 


(2.7) 


( 2 . 8 ) 


Coefficients of this expansion can be found from the continuity of p ^ an d the corresponding derivatives 
dp s ldA s and d 2 p s ld A 3 at A s = 0. After setting A s = 1, the expression ' 2.2 yields the relation between the 
pressure P and the excess chemical potential pf 1 of a fluid: 




ex 

1 


ju‘j ) ) = -ln(l-77i/0 o ) + ^ 


771/00 

\~T]il(pO 


] E (??i/ 0 o ) 2 
(l-t?l/ 0 o ) 2 


pP_rp 

0 Pi ’ 


(2.9) 


where the coefficients A and B determine the porous medium structure, and for a HS fluid in a HS matrix, 
they are as follows: 


3?7 0 t(t + 4) 9t?qT 2 

1 -po + (l-po) 2 ’ 
™ \2 


( 2 . 10 ) 


Using the Gibbs-Duhem equation, which relates the pressure of a fluid with its total chemical potential 
pi - ln(A 3 pi) + p* x as 




( 2 . 11 ) 


one derives the fluid compressibility as 


P 


( dP 

Upi 


1 


T (1-171/0) 
+ (A + 2B) 


+ (1 + A) 


t7i/0o 


(l-77i/0)(l-?7i/0o) 

(t7i/ 0 o ) 2 


+ 2B- 


(l-77i/0)(l-?7i/0o) 2 

[mifpo ) 3 


(l-J7i/0)(l-/7i/0o) 


3 ' 


( 2 . 12 ) 


After dividing the expression \2.12) by p\ and subsequently integrating it over p\, one obtains the chem¬ 
ical potential: 


Pp\ PT2 = ln(A 3 pi) - ln( 0 ) - ln(l - 771 / 0 ) + (A+ 1 ) 


0 


-In 


I- 771/0 


+ (A + 2B)- 


0 I 171 / 0o 


0 


+ 2 B- 


0 


0-00 V 1 -771 /00 

1 (p 1 / 0 o ) 2 


-In 


0-0o 1 771 / 00 

1 - 771/0 ) 


0-00 


0 2 


rln 


2 (1-77 1 /0 O ) 2 
l-7?i/0 


0-0o 1 77 r / 0o 7 

0 Pi/00 

0 - 0 o 1 — 77 i / 00 


(0-0o) 2 1 — 771 / 00 


(2.13) 


It is worth noting that the second term - ln(0) in 1 2.13 follows from the relation' 2.5 and th e corr espond- 
ing substitution /3p® - - In(0). Similarly, integration of the right-hand side of expression 1 2.12 > over p\ 
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leads to the pressure 


PP\ 
Pi 1 


SPT 2 


= - —In- 1 " 


Pi 

+ (A + 2B) 


^ +a + A] * 


<P 


1 — t?i / <po 

<P 


1 


P\ <p-<pa 

<P </> 


-In 


+ 2B- 


<P 


(p-(po 

<P <P 2 


(p-(po \l-pi/(po 

1 Vi/(po 

2 (l-77i/0o) 2 

l-i7i l( P 


Pi <P~<Po 

2(p-(p 0 


I — 771/0 
1 — 771 /0 O 
I- 771/0 1 


-In 


1 — 771 /0 O 7 


1 


(p-(po 1 — 77 i/ 0 o 


-ln- 


(2.14) 


771 (0-0 o ) 2 I — 771 /0 O . 

The expressions (2.13} and (2.14) are considered as the expression derived within the framework of the 
SPT2 approach (20). A simple analysis of (2.13} and (2.14} shows that they have two divergences at 771 = 
<p and pi = (pQ. Since <p < <po, the divergence at pi - <p occurs at lower densities. However, from the 
geometrical point of view, this divergence should appear at higher densities near the maximum value 
of the fluid packing fraction available for a fluid in a given matrix. Different corrections improving the 
SPT2 approach were proposed in (2&H22j . Here, we consider two of them, which provide rather accurate 
results in comparison with computer simulations. The first of them known as SPT2b can be derived if c p 
is replaced by cpo everywhere in ( 2 . 12 } except for the first term. In this case, the chemical potential and 
the pressure of a confined fluid are as follows: 

Pilipo 


/3yuf PT2b = ln(Ajpi) - ln(0) - ln(l - 771 /(/>) + (1 + A) 


1 — 771 /0 O 


1 (771 /0o) 

+ ~{A + 2B) — Y , 

2 (1 -77i/0o) 2 


, 2 c (mUPof 
3 (1 - Pil<po) 3 


(2.IS) 


PP\ 
Pi 1 


SPT 2 b 


Pi 


(i- —' 

I + — lnl 

1 -— | 

l <p > 

1 p 1 I 

1 0o / 


1 


1 — 771 /0 O 


A 

+ —- 


Pilipo 


2 B (pilipof 


(2.16) 


2 (1 - 771 /0 O ) 2 3 (1 — 771 /0 O ) 3 

The second approximation referred to as SPT2bl can be derived from SPT2b by removing the divergence 
at pi - (p by an expansion of the logarithmic term in (2.15} 

77(00 — 0) 


-ln(l-77i/0)s=-ln(l - 771 / 0 O ) + , - 

000(1 — 77l/0 O ) 

As a consequence, one obtains the following expressions within the SPT2bl approximation: 

,SPT 2 bl _ w . 3 „ ! ,„ fn „ ,. t , , ,, , Pll<Po 


(2.17) 


/3/ij = ln(AjPi) - In(</>) - ln(l - p x /0 O ) + (1 + A) 

, Vltyo-W (T7i/0o) 2 

000(1 — 771 /0o) 2 (1 — 771 /0o) 2 


1 — 77i/0o 
] 2 e (/?i/0o) 3 
3 (1 - 771 / 0o) 3 


(2.18) 


PP) 

Pi 1 


SPT 2 bl 


1 


(p 0 


A 

+ — 


1 - 771 / 00 (p 

Pllipo 


l*-l] 

— In [ 1 - 

—) 

l (p 1 

Pi 1 

‘Pal 


2 B (771 UPof 


(2.19) 


2(1-7 7l /0 o ) 2 3 (l-7 7l /0 o ) 3 ' 

2.2. BH perturbation theory for simple fluid in random porous medium 

The next step of theoretical treatment is connected with a consideration of an attractive part of inter¬ 
action. We consider a simple fluid with an intermolecular interaction in the form 


i'll (r): 


00 , 

un(r), 


r < ay, 
r > a 1 , 


( 2 . 20 ) 
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where u n (r) 0 is a pure attractive part of interaction. 

In order to take into account the attractive part of interaction, in this subsection we generalize the BH 
perturbation theory for the case of a fluid in a random HS matrix. To this end, we use the replica trick (Til 
according to which a system of a fluid in a matrix of unmovable (frozen) particles can be replaced by 
an equilibrium mixture consisting of the movable (annealed) matrix particles and s identical copies (or 
replicas) of a fluid. The condition is also set that the fluid replicas from different copies do not interact 
with each other, but they interact with the matrix. Such a system can be described in a standard way 
using the liquid state theories, and the properties of a fluid can be obtained by considering the limit 
s —► 0. Therefore, the Helmholtz free energy of a fluid in a matrix can be presented as (30): 

d 

F — lim —-F(s), (2.21) 

s—o ds 


where F(s ) is the free energy of the (s+ l)-component equilibrium mixture. 

Within the framework of the Zwanzig high-temperature perturbation theory 0, the first term of the 
free energy expansion corresponds to the high-temperature approximation: 

/j(F ~^ 0) -= \plP J dfg^(r)Mn(r), (2.22) 

which is of the same form as in the optimized cluster expansions (30). Fq and g^ s (r) are the free energy 
and the pair distribution function of a HS fluid in a HS matrix, respectively. 

The second correction term of Zwanzig expansion involves the three- and four-body distribution func¬ 
tions, for which it is difficult to find simple satisfactory approximations. Therefore, instead of this, we 
follow Barker and Henderson (Hi) and we write the free energy in the form: 


P(F-F 0 ) Bli 

V 


\At >/drgff 


(r)un(r) 



'd[pign S (r)] 

, dpi 



(2.23) 


where we use the same semimacroscopic arguments as in (3J0) and neglect the cross-correlation terms 
between fluids from different replicas. 

Differentiating the expressions (2.22) and (2.23) with respect to the fluid density, one derives the ex¬ 
pression for the chemical potential of a fluid: 


Ppi = Pp f+ M ,1a + /K m > 


,HTA 


,BH 


(2.24) 


where p^ s is the HS contribution of the reference system, which can be given by equation 1 2.15 > within 
the SPT2b approximation or by the equation (2.18) within the SPT2bl approximation. The first term of 
(2.23) corresponds to the HTA approximation, and it has the following form: 


PP™ = 271P Pi 


2/(pi) + pi-—J(pi) 
dp i 


(2.25) 


The contribution coming from the BH approximation as the second term of (2.23) is as follows: 

fd Pl \ ns rr ^ d (d Pl ' 


PpT = -nfipi 


dpi 


2JiPl) 'dP 


d (dp i' HS 
+ Pi/(Pi)t— 
t d Pi V dP 


dP 


HS 


+ 4 ?1 htt + p 1 


dp i 


d/(pi) d 
dpi dpi 


dpi 

dP 


T 

HS 


+pi 


f apl ) 

HS d 2 J{pi) 


t dpj 


(2.26) 


The expressions for the first and second derivatives of the isothermal compressibility wi th res pect to the 

(do p, (fin ^HS 

fluid density for the reference system, and ^ , are found using equations' 2.16 1 and 1 2.19 
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within the framework of the SPT2b or SPT2bl approximations, correspondingly, and they are presented 
in Appendix. The functions /(pi) and /(pi) are the integrals from the first and second terms of {2.23} 

OO 

Hpi) = J gu S (r)u n (r) r 2 dr, 
o 

OO 

J(p i) = J g^(r)u 2 n (r)r 2 dr. (2.27) 

o 


These integrals contain the pair distribution function g^ s (r), which is unknown for this moment and will 
be considered separately in the next subsection. 

The pressure can be calculated by differentiating the expression )2.23) with respect to the volume of 
a system or from the general thermodynamical relation: 


Pp=Pp m~Py- 


(2.28) 


The general form of the pressure is as follows: 

PP = /3P n * + p P nlfl + PP D 


(2.29) 


where P HS is the HS contribution given by ' 2.16} within the framework of SPT2b approximation or by 
[2.19) within the framework of SPT2bl approximation. The HTA term for the pressure is as follows: 


PP 


HTA 


Pi 

The contribution of BH term is as follows: 

ppBH 


■ 2nPpi 


Hpi) + pi-z—Kp i) 

dpi 


(2.30) 


Pi 


■ P Pi 


dpi 

dP 


HS 


/ dpi ) HS d n< 

+ 3p ^) T ^ ) + p P 


+P 2 i 


, 3 ( 

dpi) HS 

dpi i 

dP ) t 

l) d 

(dpi) HS 

dpi 

idpj r 


[ 3Pl ] 

HS d 2 /(pi) 

i dP 

T dp 2 


(2.31) 


2.3. The replica Ornstein-Zernike equations 


The pair distribution function g^ s (r) needed to calculate the integrals '2.27 can be obtained from 
a solution of the so-called replica Ornstein-Zernike (ROZ) equations, which were derived by Given and 
Stell [11) using the replica trick: 


hoo = coo + Po c oo ® hoo » 

h\o = cio + pocio ® hoo + pic c ® bio, 

hn = cn + p 0 cio ® hoi + pic c ® h n +PiCb® h c , 

h c = c c +pic c ® h c , (2.32) 


where the symbol ® denotes a convolution. The pair and direct fluid-fluid correlation functions are sep¬ 
arated into the connected and blocked parts 

/7n(r) = /t c (r) + h b (r), 

Cn(r) = c c (r) + c b (r). (2.33) 

As usual in the liquid state theory 116 ITTl . the ROZ equations need additional closure relations. The 
Percus-Yevick (PY) approximation is used for a HS fluid in a HS matrix considered in this paper as the 
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Figure 1. (Color online) Fluid-fluid pair distribution functions g^ s (r) of a HS fluid in a HS matrix for 
different parameters. Left-hand panel: the different matrix densities pg = porf and the size ratios of 
fluid and matrix particles t = cri /op, but the fixed fluid density pj = picr^ = 0.5. Right-hand panel: the 
different fluid densities Pp but the fixed matrix density pg = 0.3. 


reference system. As it was shown in [33; l34l . this approximation gives results for the pair distribution 
functions which are in good agreement with computer simulations. In the PY approximation for a HS fluid 
in a HS matrix, the blocking direct correlation function is zero ct>(r) = 0, thus c c (r) = c n (rj. Therefore, 
the closure conditions in our case are as follows: 

h 00 (r) = -l if r<a 0 , c 00 (r) = 0 if r>a 0 , 

h\a{r) = -1 if r<a 0 i, ci 0 (r) = 0 if r > ct 0 i , (2.34) 

hn[r) = -l if r<o\, cn(r)=:0 if r>o i, 

where opi = o + cri). 

The set of equations (2.32) in combination with the closure relations |2.34| is solved numerically using 
the hybrid Newton-Raphson procedure (35). Some of the results for gp (r) - 1 + h u (r) are presented in 
figure [l] As one can see, the function gp S (r) shows a typical behaviour for a HS fluid. In the left-hand 
panel of figure[T|it is shown that for the fixed fluid density, an increase of the matrix density, i.e., lowering 
the matrix porosity, leads to a short-range order increase. The same effect is observed if the fluid density is 
increased, but the matrix density is kept constant (figure]!] right-hand panel). There is also a comparison 
of the ROZ results with the Monte-Carlo simulations for the low and high fluid densities of a HS fluid in 
a HS matrix. It is clearly seen that they fit very well, except the contact value gn((7j) for the dense fluid, 
which is somewhat lower in the case of ROZ equations. Besides, there are two points of the contact value 
obtained from the SPT for a fluid in a matrix (32), which are a bit higher than the values of ROZ, thus 
closer to the corresponding simulation results. In this paper, for a comparison we also consider the mean 
spherical approximation (MSA) for the studied model of a fluid with the pair interaction between particles 
in the form (2.20) . Similarly to the HTA and BH approximations, to calculate the chemical potential and 
the pressure within the MSA approximation, the correlation functions are needed. For this purpose, we 
can also use the ROZ equations in combination with the corresponding closure relations: 

h 00 {r) = -1 if r <a 0 , c 00 (r) = 0 if r > cr 0 , 

hio(r) = — 1 if r<a 0 i, c w (r) = 0 if r>cr 01 , (2.35) 

h n (r) = -l if r < ai, c n (r) = c c [r) = ~/3u n (r) if r>o\. 

Again, using the replica procedure Qj]| and the general expressions of Hoye and Stell (36) for thermo¬ 
dynamic properties of an equilibrium mixture of (5+ 1) components in the limit s —► 0, one obtains the 
chemical potential and pressure of a fluid within the framework of the MSA approximation: 

Pin = Pn f+ M* SA 

P P = PP HS +/3P msa , (2.36) 
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where “HS” denotes the reference system, i.e., the HS system. For the perturbation part within the MSA 
approximation, the expressions for the chemical potential is (30): 

OO 

/3p^ SA = 2nfipi J uu(r)g n (r)r 2 dr 
0 

OO OO 

— 2.71 p\ J [cnfrJ-cftV)] r 2 dr-2npo j [c 0 i(r) -c^(r)] r 2 dr. (2.37) 

o o 

In the same way the pressure of a fluid confined in a matrix can be found: 

pp MSA l Pl = ~n Pl a\ {gn(o-J-)- [gnK)] HS } 

+ ^7rpo^ 0 {gm^o)- [^oK 0 )] hs } 

OO 

2 C du\\{r) q 

+ -^pi J gn(r)——— r 3 dr, (2.38) 

o 

where c^ s (r) and c^ 0 s (r) are the fluid-fluid and fluid-matrix direct correlation functions for a HS fluid in 
aHS matrix, gn(fj|) and gio(crf 0 ) are the contact values of pair distribution functions gn(r) = 1 + h n (r) 
and gio(r) = l + / 2 i 0 (r). 


2.4. Some calculation details 


In the present study we consider a fluid of spherical molecules confined in a HS matrix. According 
to (2.20) , the fluid-fluid interaction is decomposed into a hard-sphere part and a Lennard-Jones (12-6) tail 
following Weeks, Chandler and Andersen |6), i.e., 


j ~£, <J\ < r < v/2cti , 

1 4e [(cxi/r) 12 - (cri/r) 6 ], r > \[2o\. 


(2.39) 


Applying a hard core potential for the repulsive potential, we avoid difficulties connected with the tem¬ 
perature dependence of the fluid diameter, which requires a special treatment in the case of soft repul¬ 
sion m. The Lennard-Jones attractive tail is truncated at r - 2.5cri in order to compare our results with 
the corresponding computer simulations (28) . 

The functions 7(pi) and /(pi) used in the generalized BH theory should be calculated from the inte¬ 
grals (2.27). The dependencies of /(pi) and /(pi) on the fluid density are illustrated in figure [2] As one 




Figure 2. (Color online) The integrals 7(pi) and 7(pi) defined by (2.27| for different matrix densities. 
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can see, these dependencies are rather smooth, that is why they can be interpolated as polynomials. It 
is noticed that polynomials of the order of 9 provide a satisfactory fitting for the given functions. Since 
finally the functions I{p\) and J{p\) are polynomials, all first and second derivatives of these functions 
used in (2.25) - {2^26) and i | 2.30} - (2~3T| are calculated analytically. 

Having the chemical potential and pressure as functions of p\ (or 771 ) at different temperatures, one 
can calculate the coexistence curves of the liquid-vapour phase transition. For this purpose, we solve a 
set of two non-linear equations which follows from the conditions of thermodynamic equilibrium: 

P\{p v v T) =p ] {p\,T), 

P(p\,T) = P{p\,T), (2.40) 

where p\ and p\ are the fluid density of vapour and liquid phases, respectively. The numerical solution of 
these equations is realized using the Newton-Raphson algorithm. Thus, the liquid-vapour phase diagrams 
of a fluid confined in HS matrices are constructed. 


3. Results and discussions 


We apply the theory presented in the previous section for the description of liquid-vapour phase co¬ 
existence of a simple fluid in a HS random porous medium. However, we first consider the bulk case, 


i.e., when pj = 0. The fluid-fluid interaction v\\ (r) is taken in the form 1 2.20 1 with the attractive poten¬ 
tial Un (r) (2.40) . In figure[3] one can see the liquid-vapour phase diagrams in coordinates T-pi for a bulk 
fluid obtained within different approximations. The simulation results obtained in (28) using the method 
of grand-canonical Monte Carlo are shown for a comparison. As it is seen, the HTA approximation gives 
a good description only for low temperature and leads to the overestimation at higher temperatures. The 
BH approximation slightly improves the diagram, but still essentially overestimates the critical tempera¬ 
ture. As expected, the MSA approximation provides the best description of the phase diagram among the 
considered ones, although the critical temperature is still higher than in the simulations. 

Considering a fluid in a matrix within the model applied in our study, one should take into account 
that there is no attractive interaction between fluid and matrix particles. Thus, the only effect of the 
matrix is to confine the fluid in the void volume formed between matrix particles. Therefore, the most 
relevant parameters, which determine fluid properties, are the matrix density pjj or the corresponding 
matrix porosity (f > 0 , and the size ratio of fluid and matrix particles t - 0 \la 0 . In figure[4] we present the 
liquid-vapour phase diagrams obtained for a fluid in matrices of different densities pg = 0 (bulk), 0.046, 
0.15 and 0.30. In this case, the fluid and matrix particles are of equal sizes, i.e., t = 1. For this purpose, 



Figure 3. (Color online) Liquid-vapour phase diagram for the bulk fluid (po = 0) with the fluid-fluid in¬ 
teraction (2.20) and (2.39) , where T* = k%Tle. Solid line corresponds to BH theory, dotted line — HTA, 
dashed line — MSA, symbols — GEMC simulation results (28). 
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Figure 4. (Color online) Liquid-vapour phase diagrams for the same fluid as in figure[3] but confined in 
a matrix of different densities pg = pofff = 0.046, 0.15 and 0.30, where T* = k B Tle. Left-hand panel: 
dashed lines correspond to the HTA approximation combined with SPT2b, solid lines — HTA combined 
with SPT2bl, symbols — GCMC results taken from (28). Right-hand panel: dashed lines correspond to the 
BH theory combined with SPT2b, solid lines — BH theory combined with SPT2bl, dotted lines — MSA 
approximation, symbols — the GCMC results taken from f28| . 


we use the theoretical approaches considered in the previous section, i.e., the HTA and BH approxima¬ 
tions in combination with the reference system obtained within the SPT2b and SPT2bl approaches. The 
computer simulation results taken from (28) are presented in figure [5] for comparison. Similarly to the 
bulk case, one can see the MSA results for the considered systems (30). It is observed that for a low matrix 
density p g = 0.046, the coexistence curves calculated using the SPT2b and SPT2bl approximations almost 
coincide. With an increase of the matrix density up to pg = 0.15, the difference between the results ob¬ 
tained with SPT2b and SPT2bl becomes more distinguishable. And finally, for the high matrix density 
p* = 0.3, the diagrams differ essentially in these approximations. Moreover, the results obtained with the 
use of the SPT2b are rather anomalous, since they are far from any other approximation. This anomaly 
can be explained by the divergence contained in the expressions for the chemical potential and pressure 
of a fluid when pi — <p (SUED, and which become important at high matrix densities. To illustrate this 
problem, the dependencies of the chemical potential of a HS fluid in a HS matrix of densities p* = 0.15 
and pg = 0.3, in comparison with the results of grand-canonical Monte Carlo are shown in figure [s] It is 
clearly seen that for the matrix density pg = 0.3 the chemical potential of a fluid becomes wrong at the 
densities pj > 0.2 and tends to infinity around 0.33. The latter value corresponds to the pi ~ (f>, where 
the divergence is expected. However, this is not the case for pg = 0.15, where the SPT2b approximation 
is close to the result of SPT2bl approximation. On the other hand, it is observed that the SPT2bl approxi¬ 
mation perfectly fits the simulation results for the both matrix densities up to the highest values of fluid 
densities. Therefore, STP2bl should be considered as the best choice for the description of a HS fluid in 
HS matrix. Consequently, hereafter we restrict ourselves only to this approximation for the reference 
system. 

All the approximations HTA, BH and MSA correctly reproduce the basic trends of the behaviour of 
liquid-vapour coexistence curves of a fluid in a matrix, i.e., a decrease of matrix porosity (or an increase 
of the matrix density) leads to a critical point shift toward lower fluid densities and lower temperatures, 
simultaneously the phase diagram becomes narrower (figure [4j. Furthermore, all of the considered ap¬ 
proximations give only one critical point. A comparison of the diagrams obtained using the approxima¬ 
tions presented in our paper with computer simulations (figure[4j shows that the inclusion of the second 
term in the BH theory essentially improves the description of coexistence curves. The computer simula¬ 
tions data are in semiquantitative agreement with the theoretical prediction based on the MSA and BH 
approximation. It is worth mentioning that in contrast to MSA, which is mostly numerical approach, the 
BH approximation is more an analytical theory. 

The phase diagrams obtained for a fluid in the matrices of different densities are presented in the 
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Figure 5. The chemical potential for a HS fluid confined in a HS matrix of densities pj = po&f = 
0.15 and 0.30, and the ratio t = <ti/cto = 1.0. A comparison between the results of SPT2b (dashed 
lines), SPT2bl (solid lines) and the grand-canonical Monte Carlo simulations performed in the present 
study (symbols). 


reduced units as a plot T/T cr versus pi/pi, C r infigure[6] In such a way, it makes possible to check whether 
the theory developed for a fluid confined in a HS matrix satisfies the law of the corresponding states. As 
one can see in the figure, the considered diagrams are rather close to each other, except the case of high 
matrix density pg = 0.30. A general trend of the phase diagrams becomes broader when the reduced 
temperature decreases. 

Now, we consider the confinement effect of a matrix on a fluid by varying the size ratio of the fluid and 
matrix particles r = o\Ioq at a fixed porosity. At the same time, we compare the theoretical results with 
computer simulations data obtained in (29) using the method of Gibbs-ensemble Monte Carlo (GEMC) for 
the conventional Lennard-Jones (LJ) potential 


Pn(r) =4e 



(3.1) 


truncated at the distance r c = 2.5tJi. Since in our study the reference system is taken as a HS system, the 



Figure 6. Liquid-vapour coexistence diagrams in terms of reduced temperature T / T cr and reduced den¬ 
sity pi/pi, C r for different matrix densities pg and for the same fluid model as in figures [3j[4] All curves 
are calculated within the framework of the BH theory with the description of the reference system within 
SPT2bl approximation. 
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Figure 7. (Color online) Liquid-vapour phase diagram for a Lennard-Jones fluid in a HS matrix for dif¬ 
ferent size ratios t = (Ti/cro, hut at fixed porosity <pQ = 0.95. Theoretical predictions are given using the 
BH theory with the description of the reference system within framework of the SPT2bl approach. Solid 
lines — reference system with a hard core size d\(T) = o\. The dashed lines — corrected results with 
d\[T) defined according to the Barker-Henderson formula (3.4). Circles — GEMC results taken from (29). 


repulsive part of LJ potential is a hard core potential. Thus, according to the form (2.20) , we use 

{ oo, r < cri, 

4 e[(( 7 i/r) 12 -(tJi/r) 6 ], 04 < r < 2.5or , (3.2) 

0, r> 2.5cti. 

In figure [7] (solid lines) there is presented a comparison of theoretical results calculated for the fluid- 
fluid potential HI with computer simulations (29) (symbols) for the LJ potential o at r = 1 and t = 3/2, 
but for the fixed matrix porosity (f> 0 = 0.95 (or the matrix packing fraction 770 = 0.05). As one can see, 
the form of the phase diagrams obtained from the BH theory is rather close to the data of simulations, 
although for the both cases of r they are notably narrower. This deviation is systematic and anticipated, 
since the simulations are performed for the conventional LJ potential, which has a soft core, while we use 
the reference system as a hard core fluid. According to the repulsive potential is as follows: 

ipnlrl = I r «7„ 

[ 0, r > (Ji. 

It is acceptable to substitute this repulsive part of the LJ potential by a hard core, but instead of the 
diameter of HS particles 04 , one should take an effective diameter d\ which is somewhat smaller than <J\ 
and depends on the fluid temperature. To take into account a correct value of the effective diameter d\ 
we use one of the successful and the most popular relations proposed by Barker and Henderson for a LJ 
fluid (3): 

or 

df H (!T)= J {l-exp[-/3<pn(r)]}dr. (3.4) 

0 

Simple calculations of d^ H {T) depending on temperature show that it does not vary too heavily. For in¬ 
stance, for the temperature T* = 1.1, the effective diameter is = 0.9711 and for T* - 0.6, the diameter 
is = 0.9815. However, it still can have a strong effect on the thermodynamics of the system, hence on 
the curves of the liquid-vapour coexistence. Therefore, we need to substitute the diameter of HS particles 
£ 7 i by rf® H (r) in every place where it is needed in the expressions used for the reference system, i.e., the 
terms containing 04 or depending on it should be modified. The corrected expressions for the chemical 
potential and the pressure for a fluid in a HS matrix are used and calculated depending on the temper¬ 
ature. Again using the BH theory in combination with the SPT2bl approximation, the phase diagrams 
of liquid-vapour transition are obtained. The results of this correction is presented as dashed lines in 
figure[7] and, as one can see, the theoretical curves coincide very well with the computer simulation data. 
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Figure 8. The dependencies of the critical density cr and critical temperature r c * r on the matrix porosity 

<pa for the same model as in figures [3jj4] 


It is worth noting that there is a quicker and more efficient way to improve the present phase di¬ 
agrams. Our preliminary calculations show that coexistence curves change negligibly along the '/'-axis, 
and most deviations take place along the pi-axis. This is mainly related to the effect of an excluded vol¬ 
ume, which depends on 771 = ndfpi/6 and is overestimated in the case of the reference system with a 
hard sphere size d\ = (J\ . Therefore, to improve our results for the conventional LJ fluid, the packing 
fraction should be replaced by 77 ™ = n[d^ n {T)] 3 pi 16, and this is equivalent to the rescaling of the fluid 
density as: 


P? H = Pi 



(3.S) 


The correction made in such a way allows us to obtain the diagrams which are practically equivalent to 
those shown in figure[7] 

Finally, in figure| 8 ]we present the critical temperature T* r and the critical density p* cr as a function 
of porosity (po for the fluid confined in a HS matrix within the framework of the model considered in 
figure[4] One can see here the effects similar to those observed in figure[4]and figure[7] i.e., with a decrease 
of porosity <po, the critical point shifts toward lower temperatures and densities. In figure 8 we do not 
present the results for the model discussed in figure[7] since the results in figure[4]and figure 7 are shown 
for a different model of the fluid. However, we should remark that the critical temperature decreases with 
an increase of t at the fixed porosity, and the change of the critical density is very small. For example, 
from the results presented in figure [7]it is found that for <po = 0.95 and r = 2/3, the critical temperature 
r c * = 1.095 and the critical density is p* cr = 0.273, while for x = 1, the critical temperature is T* r - 1.031 
and the critical density is p * cr = 0.275. If we estimate this in the limit r —- 0, one can see that the critical 
temperature r c * shifts to T c * r tbulk) = 1.209 and the critical density to p\ cr - 

(tmik) C0rreS p 0n( j t 0 the values of the critical parameters of a bulk fluid. 


pHbuUO and % 


Picr^ 00 = 0.305x0.95 = 0.29. 


4. Conclusions 

In this paper, the Barker-Henderson (BH) perturbation theory is generalized for the a Lennard-Jones 
fluid confined in a random porous matrix. As the reference system, a hard sphere fluid in a hard sphere 
matrix is chosen. To describe the reference system, the extension of the scaled particle theory (SPT) is 
used, and two corresponding approximations are tested. It is shown that the SPT2bl approximation, 
which was developed recently, makes it possible to achieve a very accurate description of the thermody¬ 
namics of confined hard sphere systems. Combining the SPT approach with the BH theory, the expres¬ 
sions for the chemical potential and the pressure of a simple fluid in a hard sphere matrix are derived. 
Based on the obtained expressions, the phase diagrams of liquid-vapour transition are calculated and 
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Figure 9. (Color online) Isothermal compressibility and the derivative ^ for a HS fluid 

in a HS matrix calculated within SPT2b and SPT2bl approximations. 


compared with other theoretical approaches such as the high-temperature approximation and the mean- 
spherical approximation. A comparison of our results with computer simulation data found in the liter¬ 
ature is made as well. Different matrix porosities as well as the size ratios of fluid and matrix particles 
are considered in the present paper. The proposed extension of the BH theory for the case of a fluid in a 
matrix provides a good qualitative agreement with the computer simulations, and in some situations it 
provides a marvelous quantitative agreement, as it is observed in the case of low matrix porosities and 
at temperatures which are not very close to the critical point. The theory correctly reproduces the basic 
effects of porous media on the liquid-vapour phase coexistence of simple fluids, i.e., with a decrease of 
porosity, the critical point shifts toward lower fluid densities and lower temperatures. It is also observed 
that for a fixed matrix porosity, but for variable sizes of matrix particles, the critical temperature in¬ 
creases if the size of matrix particles becomes larger and moves to the value of critical temperature of a 
bulk fluid, while the critical density changes weakly, and in the limit r -*■ 0, it moves to the bulk critical 
value normalized by the porosity 0o . 

The approach developed in this paper can be extended to the case of more complex fluid systems in 
a confinement. In future, we plan to generalize the BH theory for anisotropic fluids in random porous 
media. Our very recent investigation 23 connected with the extension of the Van der Waals theory to 
the case of anisotropic fluids in random porous matrices shows that due to anisotropic interactions, the 
orientational order causes a competition between isotropic and anisotropic interactions, and the effect 
of a matrix can essentially modify the liquid-vapour phase diagram. 


A. Isothermal compressibility for a HS fluid in HS matrix 


Here, we present the expression for the isothermal compressibility and the derivative 

flp7 ( 7 ^),,. ^ or a fluid a matrix obtained within the SPT2b and SPT2bl approximations. Using 
simple differentiations of the expressions (2.16) and (2.19), one obtains 


dpi ] 

dP J 


SPT2b 


(p 0 l(p-\ 0 O 


+ A- 


Pi/0o 


(t?i/0o) 


(l-rji/0 o ) (l — t?! / 0o) 


+ 2 5 
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[dpiV 

dpi 

Up J 


SPT2b 


l-77l/0O 0 (l-77i/0 O ) 2 

0o /0 — 1 | 2 1 + 2?7 i/0o 

0o(1-77i/0o) 2 + 0(1 t)l / 00 ) 3 + 0o(l-?7i/0o) 4 

00 1 , . t?l / 00 


-1 
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0o _ i 
0 


+ A- 


1 — t?t /00 0 (1 — t?l/0o) 2 (1 -771/00^ 
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0o(l-77i/0o) 5 
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/dpi) SPT2bl = 1 : Vil(po t A ml $o | ^ (77i/0o) 2 

' dP 1 — 771/0 (1-77!/ 0 o ) 2 (1 — 771 / 0 o ) 3 (l-r?i/ 0 o ) 4 


_d_ (^Pi] SPT2bl _ __ 1 J 1 + 771 /0o ^ 1 + 2771 / 0 o | ^^ 77i/0o + (7?i/(/)o) 2 

apilapjj. (p{l-rhl(Po) 2 + 0o(l-pi/(/)o) 3 + </>o(l-pi/0o) 4 + </) O (l-7?i/0o) 5 

x [ 1 , *7i/0o , ^ ml<Po i2B (??i/0o) 2 |~ 2 

L1 — 77i/0 (1 — 771 /0o) 2 (1 — 771 /0 O ) 3 (1 - 771 / 0 O ) 4 


where A and B are given in equation {2.10) . 

The dependence of | and ^ Rjp-j on the fluid packing fraction 771 calculated within the SPT2b 
and SPT2bl approximations for different matrix densities pj is shown in figure [ 9 ] One can see that for 
low densities of pj, the results in the both approximations are nearly identical, but for pj = 0.30, the 
results within SPT2b show an odd behaviour for 771 larger than 0.15. 
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li^o Taxe piflMHa b HeBnopsiAKOBaHOMy nopi/icTOMy 
cepeflOBMLni: Teopm 36ypeHb Bapxepa-reHflepcoHa 

M.0. T orobko, T.M. naL 4 araH, B.I. UlMomnoxa 

iHCTMTyT cf)i3HKM KOHfleHCOBaHHX CUCTeM HAH YKpaiHM, By n. I. CBEHLtiqbKOrO, 1,7901 1 JlbBiB, YKpaiHa 

3acTOcoBaHO Teopira 36ypeHb BapKepa-reHgepcoHa (Bl") g/in BUBuem-w n/inHy /leHapfla-flxoHca b HeBnopng- 
KOBaHiR nopucTiii MaTpwqi, ctjjopMOBaHiii TBepgnMM ctfiepnuHnMn uacTHHKaMn. 3 Menuo onucy CMCTeMU Big- 
/lixy, ni<a Heo6xiflHa g/in Teopii 36ypeHb : 6y/io BnxopncTaHO po3BWHeHHS Teopii Maaina6HOi uacTnHKn (TMH). 
OcTaHHi gocsmem-ifl y po3BHTKy TMM g/w TBepgOKy/ibKOBOro nmtHy b TBepgOKy/ibKOBiii MaTpuqi g03BO/imoTb 
OTpi/lMyBaTH TepMOflUHaMitHi B/iaCTMBOCTi B TaKiii CnCTeMi i3 BMCOKOK) TOUHiCTK). TaKI/IM HUHOM, HaMl/l noeg- 
HaHO Teopifo BT 3 Teopieio TMM Ta BWBegeHO Bnpa3M pp in xiMitHoro noTeHqia/iy i Ti/icxy n/ntHy b MaTpnqi. 
BnKopncTOByKwn OTpnMaHi Bnpa3n Ta yMOBH (|>a30BOi piBHOBarn, no6ygoBaHO 4>a30Bi fliarpat/ut ra3-pignHa 
n/inHy/leHapfla-/)>i<OHCa b TBepgOKy/ibKOBiti MaTpwpi. /(oc/iiflxeHO et|)eKT nopnETOCTi MaTpuqi i po3Mipy MaTpn- 
hhi/ix uacTnHOK. rioKa3aHO : mo 3MeHLueHHfl nopucTOCTi MaTpnqi nom/ixye 3HaueHHS kphtwhho'i Tewinepatypn 
i KpnTnuHoi rycTHHM n/ntHy, pa30M 3 tum, c()a30Ba giarpaMa 3ByxyETben. Taicox cnocTepexeHO, m° 36i/ibLueH- 
Hn po3Mipy MaTpuuHnx nacTHHOK npn3B0gwTb flo pocTy KpuiTUMHOi TeMnepaTypn. 3ayBaxeHO, m 0 pe3y/ibTa™ 
Teopii y3rogxyK)Tbcn i3 gamtMH KOMn'KiTepHoro MOfle/noBam-m. 3anponoHOBaHwii TeopeTwm-H'iii nigxifl Moxe 
6yTH po3BHHyTnit go orntcy aHi30Tponm/ix pifli/iH y TBepflOKy/ibKOBiit MaTpuqi. 

KniOHOBi cnoBa: nmu-m b HeBnopnflKOBaHMX nopucri/ix cepeflOBimiax, reopin 36ypeHHH 
Bapxepa-TeHfiepcoHa, cf>a30BMpi nepexifl ra3-piflMHa, reopin MacujTa6HO'i yacmHKM 
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